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DUDLE1  .      •  w 

^VAI  =  SCHOO 

EFFECTS  OF  CAVITATION  ON  UNDERWATER  SHOCK 

LOADING  -  PLANE  PROBLEM,  PART  I 
1.   Introduction 

In  Ref.  1  possible  problem  formulations  for  predicting 
the  effects  of  cavitation  on  underwater  shock  propagation  were 
examined.   It  was  concluded  that  using  either  the  displacement 
potential  (a  scalar)  or  the  particle  displacement  (a  vector) 
as  the  dependent  variable  could  produce  useful  results. 
Ref.  2  gave  a  summary  account  of  an  unsuccessful  effort  to 
use  ADINA  for  applications  with  axisymmetric  geometry. 

The  present  report  is  a  summary  of  progress  to  date  in 
developing  a  capability  for  modeling  a  two-dimensional  fluid 
region  in  which  a  structure  is  submerged.   Examination  of 
available  codes  disclosed  no  candidate  allowing  implementation 
of  the  displacement  potential  formulation,  with  cavitation, 
for  the  fluid  and  simultaneously  permitting  an  appropriately 
coupled  structural  model.   Accordingly,  an  ad  hoc  code  for 
the  fluid  has  been  developed.   A  compatible  structural  code 
is  being  constructed.   All  codes  referred  to  in  this  report 
are  written  in  FORTRAN  IV. 

2.   Program  Development 

2. 1   Mesh  Generator  Q1GNEW) 

For  a  plane  problem  with  a  cylindric  structure  whose 
length  is  normal  to  the  plane,  the  fluid  region  is  the  portion 
of  the  plane  exterior  to  a  circle  (radius  A) .   Taking  advantage 
of  symmetry  and  truncating  gives  the  rectangular  region  of 


Fig.  1.  The  mesh  is  generated  from  a  set  of  super-elements 
as  shown  in  Fig.  2.  The  transitional  super-elements  (those 
adjacent  to  the  structure)  each  have  two  curved  sides.  The 
remaining  elements  are  square.  Note  that,  if  each  quadrant 
of  the  structural  "hole"  is  treated  as  an  element,  the  region 
is  topologically  equivalent  to  MU  layers,  each  containing 
ML  +  MR  squares.  In  the  region  of  Fig.  2  the  entry  face  is 
at  x  =  MR*(0.9*A)  and  the  exit  face  at  x  =  ML*(0.9*A).  The 
top  face  is  at  y  =  MU*(0.9*A). 

In  the  finite  element  grid  each  super-element  is  sub- 
divided N*N.   The  final  elements  are  four-noded  quadrilaterals 
using  linear  shape  functions.   For  the  transitional  elements 
the  additional  nodes  are  found  by  using  quadratic  isoparametric 
shape  functions  for  mapping.*   The  square  grid  on  which  the 
element  mesh  is  based  includes  N*(2*N+1)  nodes  which  are 
inside  the  structure.   It  is  advantageous  to  retain  these  nodes 
in  the  numbering  scheme. 

All  of  the  nodal  locations  within  the  transition  region 
can  be  derived  from  those  in  the  sector  from  0   to  45   by 
successive  reflections.   It  follows  that  there  are  only 
(3*N+l)*N/2  distinct  transition  elements  for  which  fluid  "mass" 
and  "stiffness"  matrices  need  to  be  calculated.   The  remaining 
elements  are  identical  squares  and  only  a  single  additional 
pair  of  fluid  matrices  is  required. 

A  program  (MGNEW)  has  been  written  to  generate  this 
element  grid.   Required  input  data  are  the  structure  radius  A 


^Mapping  with  isoparametric  shape  functions  is  described  in 
Ref.  3. 
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and  the  mesh  descriptors  ML,  MR,  MU,  and  N.   Output  data  in- 
clude nodal  coordinates,  numbers  of  nodes  and  elements,  and 
global  node  numbers  for  each  element. 
2. 2   Fluid  Program  (DPLPOT) 

A  program  for  tracking  shock  wave  propagation,  including 
cavitation,  has  been  developed.   Equations  for  the  displacement 
potential  formulation  used  are  given  in  Ref.  1. 

2.21  Initial  conditions 

A  subroutine  of  the  program  determines  initial  values 
of  the  displacement  potential  and  its  first  temporal  derivative . 
These  correspond  to  a  uniform  hydrostatic  pressure  plus  a 
plane  shock  wave  with  step  rise  and  exponential  decay.   The 
shock  wave  is  travelling  in  the  negative  x  direction  with 
front  at  x  =  A  at  time  t  =  0. 

2.22  Boundary  conditions 

The  governing  equation  requires  that  values  of  the 
normal  derivative  of  the  displacement  function  be  specified  at 
boundary  nodes.   Subroutines  are  provided  to  impose  a  radiation 
boundary  condition  on  the  entry  and  exit  faces.   At  the  entry 
face  explicit  provision  is  made  for  the  continuing  entry  of 
the  incident  shock  wave.   The  analytic  basis  for  these  boundary 
conditions  is  given  in  Ref.  1.   A  radiation  boundary  condition 
is  also  imposed  at  the  top  face,  again  including  explicit 
provision  for  the  passage  of  the  incident  shock  wave.   In  the 
absence  of  specific  intervention  the  solution  scheme  enforces 
a  zero  value  of  the  normal  derivative  at  boundary  points . 
Accordingly,  no  boundary  inputs  are  required  for  the  plane 
of  symmetry  y  =  0. 


2.23  Temporary  structure 

Initial  testing  of  the  fluid  program  DPLPOT  has  been 
done  with  a  simplified  boundary  condition  at  the  fluid- 
structure  interface.   This  condition  is  that  the  boundary 
pressure  remains  equal  to  the  initial  hydrostatic  pressure. 
Physically  this  corresponds  to  substituting  an  inflated 
cylindric  bag  for  the  structure.   Such  a  "hole"  in  the  fluid 
can  readily  be  made  to  induce  cavitation. 

2.24  Coefficient  matrices 

Coefficient  matrices  for  the  fluid  are  calculated  at 
the  element  level.   The  fluid  "mass"  matrix  is  lumped  (diagonal) 
rather  than  consistent.   Total  storage  requirements  for  element 
matrices  are  quite  modest  because  of  the  limited  number  of 
different  elements.   The  diagonal  system  mass  matrix  is 
assembled  in  a  single  vector  whose  length  is  equal  to  the 
number  of  nodes.   No  system  "stiffness"  matrix  is  assembled. 

2.25  Time  integration 

Numerical  integration  in  time  is  based  on  central 
difference  formulas.   The  calculation  is  organized  in  such  a 
way  that  the  nonlinearity  introduced  by  cavitation  does  not 
affect  the  fluid  mass  or  stiffness  matrices.   Because  the 
integration  algorithm  is  explicit  and  the  mass  matrix  is 
diagonal,  no  matrix  decomposition  is  required.   As  a  result, 
the  central  processor  time  is  directly  proportional  to  the 
number  of  fluid  nodes.   Processing  time  is  a  linear  function 
of  the  number  of  time  steps.   "Overhead"  time  requirements 
for  mesh  generation,  coefficient  matrix  calculation  and 
establishment  of  initial  values  are  approximately  equivalent 


to  9  time  steps. 

2. 3  Fluid  Program  (VELPOT) 

The  initial  form  of  the  program  DPLPOT,  described  above, 
required  a  computationally  expensive  damping  calculation. 
To  avoid  this,  the  alternative  of  using  the  velocity  potential 
as  the  dependent  variable  was  explored.   Pertinent  equations 
are  given  in  Ref.  1.   For  this  purpose  it  was  necessary  to 
revise  four  of  the  subroutines  of  DPLPOT.   The  alternate 
version,  called  VELPOT,  reduced  the  central  processor  time  by 
35%.   Although  VELPOT  results  closely  resembled  those  from 
DPLPOT,  they  were  not  identical.   Efforts  to  account  for  the 
discrepancies  led  to  the  development  of  two  one-dimensional 
programs,  DIS  and  VEE,  which  were  made  to  yield  identical 
results . 

Subsequent  development  of  DPLPOT  replaced  the  inefficien- 
damping  calculation.  This  improvement  eliminated  the  computa- 
tional advantage  of  VELPOT  and  VELPOT  was  abandoned. 

2.4  Structural  Program  (STRUK2) 

A  structural  program  has  been  under  development  by 
Jack  T.  Waller,  a  thesis  student  of  mine.   This  program 
utilizes  trigonometric  series  representations  of  the  loading 
and  displacements  of  an  orthotropic  shell  in  a  state  of  plane 
strain.   It  is  now  being  tested  separately,  but  has  not  yet 
been  coupled  with  the  fluid. 

3.   Results 

Results  obtained  thus  far  are  relevant  to  determining 
when  cavitation  will  occur,  how  big  the  cavity  will  become, 


and  how  long  it  will  persist.   For  this  purpose,  various 
combinations  of  shock  pressure,  hydrostatic  pressure,  and 
shock  decay  length  have  been  tested.   Effects  of  region  size 
and  shape,  element  size,  length  of  time  step,  damping  parameter, 
and  total  integration  time  have  also  been  examined. 
5 . 1   Line  Printer  Plot 

The  program  provides  for  both  tabular  and  graphical 
output,  but  the  latter  has  been  the  more  useful  in  this  phase 
of  the  investigation.   It  is  a  line  printer  plot  of  fluid 
nodal  pressures,  repeated  at  selected  time  intervals.   For 
the  purpose  of  this  plot  the  nodes  of  the  structural  and 
transition  regions  are  mapped  into  the  topologically  equivalent 
rectangular  array. 

Fig.  3  shows  the  arrangement  of  the  line  printer  plot. 
The  left-hand  column  is  the  plane  of  symmetry  (x  axis  of  Fig.  1) 
and  the  right-hand  column  is  the  "top"  face.   The  entry  face 
is  the  first  row  and  the  exit  face  is  the  last  row.   Along  the 
left-hand  edge  is  a  rectangle  filled  with  X's,  4  columns  wide 
and  7  rows  high.   The  X's  denote  dummy  nodes  inside  the 
structure.   The  rows  immediately  above  and  below,  and  the 
column  to  the  right  (all  marked  H)  contain  the  17  structural 
nodes . 

The  pressure  map  of  Fig.  3  is  for  hydrostatic  pressure 
0.2  MPa  and  shock  pressure  2.0  Mpa  with  a  decay  length  of  15  m. 
The  pressure  ranges  for  the  mapping  characters  are  shown  at 
the  right.   Since  the  map  is  for  t  =  0  there  are  no  pressures 
below  hvdrostatic. 
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Fig.    3.      Initial    pressure  map 


3.2   Cavitation  Example 

Fig.  4  exhibits  a  sequence  of  pressure  maps  at  16  ms 
intervals  based  on  the  initial  conditions  of  Fig.  3.   Note 
that  the  cavity  develops  rapidly  in  the  first  two  frames, 
shrinks  in  the  next  two,  and  is  completely  closed  by  30  ms . 
The  succeeding  frames  in  which  the  cavity  opens  again  represent 
the  first  instance  in  which  this  behavior  was  encountered. 
It  has  been  established  by  a  longer  duration  run  that  the 
cavity  does  not  open  again  after  144  ms . 
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